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Abstract: This paper presents an extension of the matrix element method to next-to- 
leading order in perturbation theory. To accomplish this we have developed a method 
to calculate next-to- leading order weights on an event-by-event basis. This allows for the 
definition of next-to-leading order likelihoods in exactly the same fashion as at leading 
order, thus extending the matrix element method to next-to- leading order. A welcome 
by-product of the method is the straightforward and efficient generation of unweighted 
next-to-leading order events. As examples of the application of our next-to-leading order 
matrix element method we consider the measurement of the mass of the Z boson and also 
the search for the Higgs boson in the four lepton channel. 
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1. Introduction 

The continued successful running of the LHC is already resulting in an impressive data set 
with which to test the Standard Model (SM). One of the main aims of the experimental 
program is to observe the mechanism behind electroweak symmetry breaking, for which the 
postulated Higgs boson is a theoretically well- motivated example. Using the 5 fb _1 data 
set the LHC has tightly constrained the mass of the Higgs boson, whilst also providing 
tantalising hints in the low mass region (~ 120 — 125 GeV) [1, 2]. Present analyses often 
use data driven techniques for background estimation with an emphasis on accurate signal 
modeling, for instance in the diphoton Higgs search [3, 4]. Whilst this is a sensible strategy 
for searches, after discovery an accurate modeling of both signal and background will be 
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required in order to confirm the exact properties of any new particle, such as its spin and 
couplings [5, 6]. In addition to Higgs searches, precision measurements in the electroweak 
sector of the SM could also provide valuable insight. By measuring top quark properties and 
electroweak gauge boson couplings, potential new physics contributions can be constrained. 
A recent example, that exhibits some tension with the SM, is provided by measurements 
of the top quark forward-backward asymmetry at the Tevatron [7, 8]. 

There are many methods available for performing studies of particle properties, for 
instance for measuring their masses or investigating their interactions. Among these, the 
matrix element method (MEM) stands out since it is sensitive to all the available kinematic 
information for each individual event. Originally pioneered at the Tevatron [9, 10], the 
MEM has proven extremely useful in the the top sector [11-16]. Recently the method 
has been used to observe single top production [17-20] and to provide evidence for top 
quark spin correlations [21]. The MEM has also been used to try to improve searches for 
the Higgs boson in the associated production channel [22]. At the LHC the MEM is also 
beginning to be used, for example in the measurement of the electroweak mixing angle at 
CMS [23]. 

The popularity of the MEM is based on its ability to utilize the theoretical prediction 
from the matrix element, retaining all the hard scattering correlations. For each experi- 
mental event, the MEM assigns a probability that it can be described by a given theoretical 
model. In this way one can produce a likelihood that the theoretical model describes a 
particular set of data. Matrix elements (at tree-level) are relatively straightforward to 
calculate and automated tools for this purpose have been available for several years [24- 
28]. Indeed, the application of automated tools to the MEM was previously considered in 
ref. [29]. However, a serious limitation of the method is that it has so far been defined only 
at leading order (LO). For the precision studies that will become possible with the wealth 
of data at the LHC, it is crucial to extend and adapt the method such that it is defined at 
higher orders. An implementation of the method at next-to-leading order (NLO), the de 
facto standard for most theoretical predictions at the LHC, is required to put the MEM 
on a solid theoretical footing and elevate the method to being a robust analysis tool. 

The absence of higher order corrections in current implementations of the MEM is 
easily understood. It is not immediately clear how to use existing NLO calculations to 
associate a NLO weight with a given exclusive experimental event. This is primarily due 
to the fact that NLO calculations include contributions from both loop and bremsstrahlung 
diagrams, that must be integrated over different physical phase spaces. As such there is 
no clear one-to-one map between an exclusive event, containing a finite number of objects 
with measured properties, and a NLO weight. Addressing this very issue is the principal 
goal of this paper. 

We therefore present a method of calculating NLO weights suitable for use with the 
MEM approach. As a welcome by-product, the method also provides a procedure for 
calculating unweighted NLO events. As a first step, in this paper we consider only the 
production of colour neutral final states. This ensures that at NLO the real phase space is 
associated with radiation from initial state partons only. We thus postpone the treatment 
of final state jets at NLO to a future publication. 
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This paper proceeds as follows. In section || we first introduce the MEM at LO and 
discuss its use in experimental analyses. Section |3| explains our extension of the MEM to 
NLO and discusses the generation of unweighted NLO events. In section Q we validate 
the code using MCFM [30-33] and Pythia [34]. Section |5] is devoted to an application 
of immediate phenomenological interest, namely the search for a Higgs boson in the ZZ* 
decay channel to four leptons. Finally in section ^ we draw our conclusions. The appendices 
describe the generation of the phase space in more detail and discuss the modifications to 
the usual dipole subtraction procedure that are required in our approach. 



2. The Matrix Element Method at Leading Order 

In this section we define the MEM at LO and discuss how it may be used in experimental 
analyses. 

2.1 Overview of the MEM 

We begin by assuming that one wishes to measure a model parameter Q, using an experi- 
mental data set {x} that contains N events x;. One method to determine the best-fit value 
of ft is to construct a probability density function in which each event is weighted by the 
LO scattering probability computed with the parameter f2. The resulting probability den- 
sity function associated with a single event x, for a given 0, can be written schematically 
as. 

P(x|fl) = 4o / dx a dx b dyJ2 h{Xa)fj{Xb) fig(p a>P 6,y) Wfry) . (2.1) 

do J ^ X a X b S 

In this equation fi(x a ) and fj{x b ) represent the parton distribution functions for partons 
of flavours i and j possessing momentum fractions x a and x b of their parent hadrons. 
BniPaiPb, y) is the LO scattering probability with partons i and j in the initial state. The 
hadron collision takes place at a centre of mass energy y/s while the flux factor entering in 



the denominator of Eq. (2.1) is the partonic centre of mass energy squared, s ab = x a x b s. 

An experimental event x is by definition a detector level event, whilst the scattering 
probability is computed theoretically at the level of partons. Therefore in order to correctly 
use the scattering probability as a probability density function one must include effects that 
model this discrepancy. The transfer function W(x, y) relates a detector level event x to a 
particle level event y that can be used to compute the scattering amplitude. This transfer 
function, dependent on the specifics of the experimental set-up, takes account of factors 
such as limitations on the energy resolution and acceptance of the detector. The transfer 
function is constructed such that it is itself a probability density function, 



/ 



cfyW(x,y) = l. (2.2) 



Finally, the factor an is the total cross section for the process for a specific choice of 
thus ensuring that the probability distribution is properly normalized to unity. 
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Once the probability density function 'P(x|J7) has been computed for each event x, it 
is straightforward to compute a likelihood for the data set as a whole. For the data set {x} 
with N events, the likelihood function £({x}|f2) for a given parameter Q is defined by, 



N 

c({ x }\n) = f(N)Hv(^\n)- (2.3) 

Here f(N) is a normalisation factor related to the overall number of events in the data 
set. In most analyses one is interested in comparing two hypotheses, either in the form of 
a likelihood ratio, or more commonly by comparing the difference of two log-likelihoods. 
Therefore in most practical applications the explicit form of f(N) is unimportant. This is 
the case for all the examples that we present here and, as such, we will simply drop the 



factor f(N) in Eq. Q. 

By construction, the value of the likelihood function will be larger for theories that 
describe the data better. The best fit corresponds to the parameter choice SI that maximises 
C (and hence also log C). In the region of the maximum - and as long as the data set is large 
enough - departures from the maximum value of the likelihood can be simply interpreted 
in terms of standard deviations from the best fit. Since we consider a single parameter 
fl, the likelihood can be described by a parabola in the region of the maximum (see e.g. 
ref [35]) and standard deviations (here represented by na) from the observed maximum 
can then be defined by, 

logC\ na = logC max -n 2 /2. (2.4) 

In our examples we will use this to define one- and two-sigma confidence levels for our 
results, although we stress that our studies do not include detector effects and are thus 
only for the sake of illustration. 

2.2 Leading order formulation 

We now return to the probability density function, Eq. ( |2.1| ). We recall that at this order, 

Bg(Pa,Pb,y) = \M% {0 \ Pa , Pb ,y)\ 2 , (2.5) 

where Mq' (y) is the leading order matrix element for the relevant process with initial 
state partons i and j. In this paper we make the simplifying assumption that the events 
completely specify the final state particles so that, for example, we do not consider events 
containing neutrinos. For a Born point p the constraint of momentum conservation fixes 
the values of the parton fractions x a and Xb- By convention we position the incoming 
particles along the z-axis in the lab frame and then use the momentum conserving delta 
function between the n final-state particles {pi, . . . ,p n }, 

n 
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Figure 1: The generation of the Born (and virtual) phase space from a given experimental event. 
The left hand side depicts a collision that results in the production of a colour neutral final state 
(represented here by four leptons in red) that do not balance in the transverse plane. The resulting 
imbalance (X, in blue) represents the remaining event which is not modelled in the Born matrix 
element. We apply a Lorentz transformation such that X has no components in the transverse 
plane, with the remaining longitudinal and energy components absorbed into the colliding partons. 





to find the relations, 

x a -x b = [ Vpf ] , x a + x b = I VE, I . (2.7) 

However, matching an experimental point p to the LO kinematics (p) is a challenge. In 
particular, any event will always contain additional radiation that is not modelled by the 
leading order (Born level) matrix element. In order to proceed we shall define a four vector 
X, that balances the momenta of the final state particles. This is illustrated schematically 
in Fig. |l] and expressed through the equations, 

n 

X = -J2Pi- (2-8) 

i=l 

The Born matrix elements, with the beam directions consistently along the z-axis, are 
only defined for X x = X y = 0, i.e. when there is no px imbalance between the final 
state particles 1 . Therefore, in order to ensure that the experimental event has a well- 
defined interpretation as a Born level phase space point we need to remove the transverse 
components of X. This can be achieved by applying a Lorentz transformation A(X) on 
the momenta p in the event to arrive at a frame in which the transverse components of X 
are zero, 

n n 

P >t = ^ v (x)pt with Y,Pi=Y,p" = °- ( 2 - 9 ) 

i=l i=l 

As desired, the phase space point p is now of the correct form to be used in a Born level 
matrix element. For a given transformation, the momentum fractions x a and x b are then 
related to the transformed momenta p through the relations in Eq. (^?j). However, we note 



1 Attempting to evaluate a LO matrix element with a phase space point that does not conserve momentum 
is ill-defined. The exact weight obtained depends on which kinematic invariants one has chosen to use in 
the expression for the matrix element. 
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that Eq. (|2.9| ) does not specify a unique transformation. We can define multiple transfor- 
mations that result in X x = X y = and that yield different longitudinal components of p. 
In other words x a and Xb are frame-dependent quantities determined by the boost choice 
and it is only the product x a x b that is Lorentz invariant. Therefore in order to produce a 
sensibly defined weight for each event we must integrate over this unobservable degree of 
freedom. 

To illustrate these ideas in more detail we begin with the usual definition of the total 
cross section for the production of n massless final state particles, 

ah° = (27T) 4 - 3 " / dx a dx b J] (^) • / ' U " )/ ' ( - r '- j Bg L + Pb - jr P )l .(2.10) 

J m=l V ltjm / XaXbS \ i= l J 

Here we have suppressed the dependence of B on the kinematics and the summation over 
i and j for clarity. We wish to factorise Eq. (2.10) into two pieces, one representing initial 
state production and the other the decay of a heavy object into the final state particles. 
To this end we define Q = p a + p b and insert the operator J dQ 2 5{x a x b s — Q 2 ) = 1, 

a{f = (2vr) 4 " 3n J dx a dx b dQ 2 5(x a x b s - Q 2 ) 



For the remainder of this paper we will define the phase space element associated with the 
final state particles as, 

= (2n)^dQ 2 f[ (^)^ (4) (q ~ £ Pm) ■ (2.12) 

m=\ \ m / y m=1 J 

Using this definition we see that, 

dx a dx b dx8(x a x b s - Q 2 )— — — B^(p a ,p b ,x.) ■ 

X a Xl)S 

= j dxCij(Q 2 ,xi,x u )B^(p a ,pb,x). (2.13) 

This separation is convenient since B l Q(p a ,p b ,x) is Lorentz invariant and need only be 
evaluated for a single phase space point. The process independent integration over boosts 
is given by, 

f fi(Xa)fj(xb) 
£ij($abi X^i) — / dx a dxb ${x a XbS Sab) 

J x a x b s 
= f t ,W°-», (2.14) 

J xi SX a S ab 

where in the second expression we have made the dependence on the upper and lower 
bounds explicit. 
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This factorisation in terms of initial and final state variables is exactly what we require 
to build our probability density function for the MEM since the experimental input is always 
a final state phase space point x. We can define Eq. (|2.1| ) more formally as, 



^( x l^) = -To [ d y^ij( s ab,xi,x u )B^(p a ,p b ,y)W(x,y). (2.15) 



For a completely inclusive description of the final state, Eqs. ( |2.14 ) and ( |2.15 ) are suf- 



ficient. However, realistic applications require transverse momentum and pseudo-rapidity 
cuts in order to define fiducial regions of the detector. It is therefore useful to consider the 
forms of the lab frame transverse momentum (p l ji b ) and pseudo-rapidity (77'°*) under the 
application of a given longitudinal boost parameterized by x a . 

The four-momenta of all the particles depend on the boost parameter - the initial 
state momenta p a (x a ), Pb(x a ) and the momentum of particle i in the final state, Pi(x a ). 
However we note that invariant masses, Sjj = 2pi(x a ) ■ Pj(x a ) cannot depend on the boost 
and may therefore be evaluated using any choice of boost parameter. The lab frame 
transverse momentum and pseudo-rapidity are defined in terms of such invariants and the 
boost parameter x a by, 



p lab,i = JSaiSib = 1 bg / X^Sib\ _ 

V s ab 2 V s ab Sai ) 

From these equations we see that p l ^ b,t does not depend on the boost parameter and 
therefore cuts on this quantity can be performed outside the boost integration, i.e. in 
Eq. ( 2.15| ). On the other hand, rj lab ' 1 depends on x a , so that cuts on the lab frame pseudo- 



rapidity should be included in Eq. (2.14). These cuts constrain the range of allowed boosts, 



i.e. the integration limits xi and x u are fixed by |r/ maa; |. 

In summary, by boosting an event to a frame in which the final state is px-balanced 
we have recovered Born kinematics and can assign a likelihood to the event uniquely. 
Frequently in the next sections we will refer to these frames, in which the Born event is 
well defined, as the "MEM frame" . As we have discussed, this definition is only unique in 
the transverse plane and the "MEM frame" is actually a set of equivalent frames connected 
by longitudinal boosts. 

For the remainder of the paper we will make a simplification by assuming a "perfect" 
detector, i.e. the transfer function is equal to W(x, y) = 6(x — y). This assumption is only 
valid for well-measured final state particles such as leptons and therefore as examples we 
only consider ZZ — > At and Z — > £ + £~~. We stress that non-trivial transfer functions do 
not pose any conceptual problems for our method and only entail additional integrations. 
Taking this simplification and the integration over the longitudinal boost into account, 



Eq. (|2.15[) becomes, 



V(x\n) = -^Cij(s ab ,x h x u )B l ^{p a ,p bl yL) . (2.17) 



The above equation defines the LO probability density function for the MEM. We recall 
that B^(p a ,pb, x) represents the Born Matrix element squared, \-M^ \p a ,Pb, x)| 2 and that 
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era represents the fiducial cross section, calculated using cuts in the lab frame. We define 
the following quantity, 

B n (x) = Cij(Sab,Xl,X u ) &g(p a ,Pb,x) , (2-18) 

and observe from Eq. ( |2.13| ) that J dxBQ(x) = (7q°. We can thus simplify Eq. ( p.!7| ) to, 

p(x|ft) = 4o^( x )- ( 2 - 19 ) 

This formalism will prove useful in the following section when we extend the MEM to NLO. 

Using the techniques outlined above we have defined a procedure that takes an observed 
final state, Q + X and relates it to a LO model for the process, p a + Pb — >• Q- Specifically, 
given an arbitrary amount of additional radiation we create a phase space point that 
recovers the Born kinematics, at the cost of introducing an integration over the longitudinal 
degree of freedom. 

Clearly this model will be better for events in which the momentum imbalance X is 
small, rather than events in which X is kinematically relevant, i.e. in the presence of one 
or more additional jets. When additional jets are present one has three options. The first 
option is to simply apply the LO model presented above, boosting the jet into the initial 
state. Since in general one expects this method to be rather sensitive to the amount of 
radiation, i.e. the transverse momentum of the jet, it is prudent to check the validity of 
this approach by also considering smaller data sets obtained by applying a jet veto. If there 
are sufficient events, restricting the data set by imposing a strict jet veto is preferred since, 
by ensuring that no additional hard jets are present, one can be confident that the LO 
model works reasonably well. We shall present an example of applying such a jet veto in 
section ^|. The second option is to use a LO calculation that already contains an additional 
jet, i.e. p a + pb — > Q + jet + X. In this case the extra radiation is well modelled but the 
MEM must be extended to include a systematic treatment of jets. In this paper we will 
not consider this option further. 

Finally, one may try to systematically improve the MEM in an attempt to model 
the additional radiation. This is the approach discussed in Ref. [36], with reference to 
initial state radiation. Instead one may incorporate such effects by extending the MEM 
to NLO. Since a NLO calculation includes the radiation of one additional parton, a first 
approximation of the effects of further radiation is made at this order. In the next section 
we will illustrate how this may be achieved within the MEM framework. 

3. The Matrix Element Method at Next-to-Leading Order 

In this section we define the MEM at NLO and, as a by-product, discuss how one may 
generate unweighted events at NLO. 

3.1 Going beyond LO: Defining NLO on an event by event basis 

The goal of this sub-section is to illustrate how to extend the MEM to NLO in perturba- 
tion theory. However this is not a simple task since in a normal NLO calculation virtual 
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and bremsstrahlung events live in separate phase spaces, their only communication being 
through a regularising subtraction scheme. Instead of following this procedure, we need to 
reorganise the calculation such that it can provide a NLO weight for a given Born event, 
with the sum over the event weights recovering the usual NLO cross section. To do this we 
begin by assuming that our event has been rendered in the MEM frame using the procedure 
described in the previous section. We note however that the procedure we will outline in 
this section is not useful solely for extending the MEM to NLO. We are creating a method 
for producing a NLO cross section from a series of Born phase space points, a procedure 
that may have broader applications than are presented here. 

Given the phase space point x = p±, . . . ,p n where the final state momenta are those 
of the identified final state particles, we can define the NLO corrections by, 

dg " W =jRn (x) + Vh(x). (3.1) 
ax 

This follows the usual separation of the NLO calculation into two pieces, each of which 
is associated with a different phase space. We stress though that here the separation has 
been performed for a fixed Born phase space point, x. The definition of the term associated 
with the virtual corrections is straightforward since it is defined in the same phase space 
as the Born contribution. Explicitly, we can define Vh(x) as, 



Vh(x) = Cij(s ah ,x h x u )\ B%(p a ,p b ,x) +V^(p a ,P6,x) 

+ ^ dz(v m (z,x.)®L m (z,s ab ,xi,x u )j B$(p a ,p b ,x). (3.2) 



Here the first term represents the combination of the Born matrix element Bq and the 
one-loop Born interference term Vq = 2Re|.A/f[^* Mq *| (where the dependence on the 
initial state partons has been suppressed). This is coupled to the same boost function, 
Cij as was defined at LO. In our approach we have followed the NLO implementation of 
MCFM and used the dipole subtraction procedure of Catani and Seymour [37] to handle 



the singularities in the virtual and real calculations. The final term in Eq. (3^) contains the 
integrated subtraction terms, T) a , introduced in this formalism. Since we are considering 
initial state singularities the integrated dipoles depend on a convolution variable z. This 
variable is convoluted with the boost function to create three structures, 



r r r /"*" a M x a/ z)fj{s ab /(sx a )) f Xu fi(x a )fj(s ab /(zsx a )) 
C = C, A = / dx a J - , £2= dx a J - , 

J xi zsx a s ab j Xl zsx a s ab 



(3.3) 



In Eq. (3.2) the sum over these convolutions is given by to. 

Using Eq. (|3.2| ) we are able to define an event by event finite weight associated with 
the Born plus virtual contributions. Our remaining task is thus to define -Rfj(x) such that 
there is no double counting of events. In other words we must ensure that the integration 
of Eq. ( |3.lD results in the total NLO cross section ((Jq LO ). One way to ensure this is to 
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use a forward branching phase space generator (FBPS) [38] to construct the real phase 
space. Starting from the Born phase space point, p a +pb — > Q the FBPS generates the real 
radiation by branching one of the initial state momenta to produce the real phase space 
point Pa+Pb Q + Pr- In the following we will use the hatted notation to indicate a Born 
phase space point, whilst the un-hatted momenta represent the real phase space point. 

The phase space generator needs to integrate out all initial state radiation within the 
constraints of fixed momenta of the identified final state particles (and, if required, the jet 
veto). We show in Appendix |A] that this can be achieved using a FBPS generator defined 
by, 

d$(p a +Pb ^Q+Pr)= d<S>{p a +p b -> Q) X d<S> FBPS (p 

aiPbiPr-) * veto 3 

(3.4) 

where 6 vcto (optionally) vetoes events that generate an additional jet. At NLO the jet veto 
cut is simply, 

0vetofrr) = 9 [p l f{p r ) < PT in (jet)] , (3.5) 

where p l ^ b {p r ) is the laboratory frame transverse momentum (calculated using Eq. ( |2.16| )). 
Note the initial state brancher is necessarily an antenna brancher since it ensures that the 
initial state partons remain massless. The form of the FBPS generator, in terms of the 
kinematic variables p a , pb and p r , is, 

d® FBPS (p a ,p b ,p r ) = — ( — J dt ar dt rb d(f> , (3.6) 
(2tt) 6 \s ab J 

where t xy = (jp x — p y ) 2 and d(j) is a rotational degree of freedom about the z-axis. The 
explicit construction of the momenta p a , p b and p r in terms of the integration variables 
is detailed in Appendix [A]. The phase space weight corrects the flux factor due to the 
resulting emission of an extra parton. 

Finally, we observe that the forward brancher must by necessity change the initial state 
momenta. This means that for bremsstrahlung events the values of p^ b will depend on the 
branching momentum p r . Thus although the four momenta of the final state particles are 
fixed in the MEM frame the value of the p l £ b observable changes dynamically. In other 
words a single event with fixed MEM frame four momenta corresponds to a range of p!f b 
values. Using the FBPS we can now explicitly define -Rn(x) as, 

Rn(x) = J d^ FWS (p a ,Pb,Pr)(^ij(Sab,Xi,X u )TZ^(p a ,p b ,X,p r ) 

-Y,£v(sab,xr,x™)D m (p a ,p b ,p r )B%(p a ,Pb,x)y (3.7) 



In the above we note that the boost integral is defined for a given branching, since each 

^(Pa,Pb,^,Pr) = 



branching generates a new s ab . The quantity TZ l A(p a ,Pb,^-,Pr) = (p a ,pb,x,p r )| 2 is the 



Born level matrix element with one additional parton. Finally, D(p a ,pb,p r )BQ(p a ,p b , x) 
represents the subtraction terms that cancel the soft and collinear divergences which occur 
when p r is unresolved. A couple of observations are in order in regards to the dipole 
pieces. We note that, since the dipoles must provide a pointwise cancellation, the boost 
function inherits the same s a b as in the real boost function. However the underlying Born 
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matrix element must be evaluated using the original Born s a b in order to have a one-to- 
one correspondence with Eq. ( |3.2| ). This also fixes the integration limits, 



Eq. (|3.7|). We discuss the exact modifications to the usual dipole subtraction scheme in 
Appendix |B|. 

We are now in a position to build our scattering probability accurate to NLO, based 
on the quantities Vh(x) and Rq(x) that we have defined in Eqs. ( |3.2[ ) and ( |3.7D above. The 
NLO probability density function associated with the event x is, 

V ( x \ n ) = ^Wlo (W) + flnfr)) • ( 3 - 8 ) 
This equation defines the MEM at NLO. 
3.2 Generating unweighted events at NLO 

A welcome by-product of the method outlined in the previous sub-section is its ability to 
generate unweighted events at NLO. In this section we outline how this is possible and in 
later sections we will use the technique to generate samples of unweighted events that can 
be used to test the MEM. 

Our starting point is Eq. (|3.l| ), in which we explicitly separated the NLO calculation 
into real and virtual contributions. We define the inclusive phase space spanned by the 
Born processes as 3>, which we can separate into two regions. Region I is the part of the 
inclusive phase space, that is populated by the LO calculation under the lab frame cuts. 
Region II is the remaining part of the inclusive phase space, in which the LO calculation 
does not contribute. 

We focus first on region I. Since the LO contribution is non-zero we can write a point 
by point K-i actor as follows, 

= do NLO fda LO \- 1 
dx \ <ix J 

" — • (3 ' 9) 

This quantity is not positive definite since one can construct phase space points for which 
Ki(x) < 0. However, these correspond to regions in which the NLO calculation is un- 
physical. More specifically, it is possible to choose a renormalisation scale such that the 
differential cross section becomes negative. Typically this occurs because the choice of 
renomalisation scale is widely separated from the typical scale of the event. In general if 
a sensible scale choice is used then Kj(x.) > 0. In order to ensure that Kj(x.) > it is 
sufficient to check that the NLO differential cross section is positive in all observables. One 
can then create weighted NLO events in this region by generating a Born phase space point 
and recording both the Born weight, Bq(x) and the -fT-factor, Kj(x.) for that point (as well 
as the phase space weight associated with x). If the calculation is completely inclusive, i.e. 
no cuts are applied and region II is empty, then an unweighted NLO sample can easily be 
obtained by unweighting the combination of Kj(x), Bq(x) and the phase space weight. 
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In region II there is no if-factor since the LO cross section is zero. In this region the 
virtual contribution, and all of the terms associated with the subtraction procedure, are 
zero since they occupy the Born phase space. Hence K/j(x) is positive definite since it 
only corresponds to the LO process with an additional parton, 



Therefore in region II we construct our weights as a combination of the phase space weight 
associated with x and Kjj(x). 

By combining regions I and II we have weights that span the entire phase space and 
which are positive (with the caveat that the total NLO differential cross section should be 
positive everywhere). Although the events all have the structure of a Born phase space 
point, the sum over the associated weights results in the NLO cross section. We stress 
that the events found in region II are those in which the Born contribution is zero due to 
fiducial cuts and not a kinematic cut off. For example if one demanded a leptonic pr cut of 
15 GeV then region II would correspond to pt < 15 GeV. On the other hand, if the lepton 
had some natural cut off (for example, px > mz/2) then this region is already excluded 
from the inclusive Born phase space, <J>. Using the weighted sample described here one can 
produce unweighted events in exactly the same fashion as one does at LO. 

4. Validation 

In this section we present a simple validation of the method outlined in the previous section, 
focussing for simplicity on the production of lepton pairs at the LHC, pp — > Z/j* — > 
£ + £~. In the first instance we study physics in the MEM frame, comparing predictions 
for observables in this frame with the more familiar ones obtained in the lab frame. For 
this exercise, we investigate parton level calculations (at LO and NLO) and Pythia [34]. 
The use of Pythia is a valuable test of our method since it contains the effects of a parton 
shower, underlying event and hadronisation in its output. After this study we present 
a simple comparison of the MEM method at LO and NLO, in the context of a Z mass 
measurement. 

4.1 Physics in the MEM frame. 

We begin by recalling the definitions of lab frame quantities (pr and rj) that we use to 
apply cuts in the MEM frame, 



In passing we note that, although it is not needed in the cases that we will discuss here, 
we can also easily define lab-frame azimuthal angle and rapidity differences in a boost- 
independent fashion, 



(3.10) 




(4.1) 




(4.2) 
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These definitions would be useful for more complicated processes that include jets or that 
require the application of an isolation procedure. We will also consider the MEM frame 
transverse momentum, which is defined in a more familiar way, 

PT m = )/(Pf ) 2 + (^) 2 , (4-3) 

where, of course, the four-vector p^ is explicitly in the MEM frame. The MEM frame has 
no unique definition of pseudo-rapidity for a given event, since there are multiple frames 
connected by longitudinal boosts. 

We now wish to study the behaviour of different quantities in the lab and MEM frames. 
We apply very loose cuts, namely we only require that the leptons lie in the invariant mass 
window, 

80 GeV < m e+i - < 100 GeV . (4.4) 

We generate LO and NLO parton level events using MCFM and more exclusive particle- 
level dilepton events using Pythia. In Fig. Q we compare the results from the lab and MEM 
frames for the quantities p l j: b , p^ E and mu- 

In Fig. H|(a) we see that, as is necessary, the invariant mass of the lepton pairs is 
identical in both frames. A more interesting quantity is the frame-dependent px of the 
positively charged lepton, £ + , shown in Fig. |2](b). At LO (parton level) the two quantities 
are the same because for pure LO results the final state has zero net transverse momentum 
and thus the MEM and lab frames are identical. As soon as this simple picture is broken 
the two frames are no longer the same and the pt distributions differ. This is apparent in 
both the showered and NLO results. For the NLO and shower predictions it is possible, 
by radiating additional particles, for a lepton to have lab-frame pt greater than mz/2. At 
LO this is not kinematically accessible, modulo small width effects. This is demonstrated 
in the lab frame pt predictions for the NLO and showered results, shown in Fig. ||(b), that 
produce a high pt tail. The MEM frame, however, requires that the event be boosted back 
to a Born topology. As such, the high px region is not present in this frame. Since the 
overall normalisation is fixed by the total cross section these events are manifested at lower 
values pt, with the region around showing a considerable enhancement relative 

to the lab frame. 

In Fig. H we directly compare the different theoretical predictions - at LO, NLO and 
using Pythia - in both the lab and MEM frames. It is clear that the predictions in 
the MEM frame are very similar with respect to each other, with both LO and Pythia 
predicting a slightly softer spectrum relative to NLO. We note that the shape differences 
between NLO and the other predictions is consistently of order 10% or less. In the lab 
frame, however, there are significant differences between the predictions, in particular in 
the region px > mz/2- From this discussion we conclude that the MEM frame possesses 
some very nice features. In particular the differences with respect to the LO prediction 
(from either shower or NLO) are consistent with naive estimates of higher order QCD 
effects, suggesting good perturbative control. The main reason for the convergence is that, 
in the MEM frame, kinematic ranges of observables are not extended beyond their LO 
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Figure 2: Comparison between the lab and MEM frame predictions from the NLO calculation of 
MCFM (left) and Pythia (right) for the process pp — > Z/j* — > £ + £~ . In (a) we plot the invariant 
mass distribution of the two leptons and in (b) we show the pr of the positively charged lepton. In 
each plot the lab frame quantity is shown in black (dashed), while the MEM frame result is in red 
(solid). 



boundaries. Since any such extension beyond the LO region is necessarily sensitive to 
further higher order corrections, the elimination of this aspect of the calculation should be 
seen as an advantage for the MEM frame. 
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Figure 3: Comparison between MCFM (LO and NLO) and Pythia in different frames. On the left 
hand side Pj, is plotted in the MEM frame, whilst on the right hand side the lab frame equivalent 
is plotted. Predictions are normalised by the total cross section (or number of events in the Pythia 
case). 



4.2 Validating the MEM: measuring mz 

In the previous sub-section we have used MCFM, representing a traditional approach to 
NLO calculations, to generate lab frame events that are then transformed into the MEM 
frame. As described in the previous section, for the the extension of the MEM method to 
NLO it is easiest to work directly in the MEM frame. We have modified MCFM accordingly 
to incorporate the phase space generator and approach described in the previous section. 
In addition to the implementation of the FBPS, the code has been constructed such that 
a NLO weight can be ascribed to an individual event in the MEM frame. 

A simple test of our implementation of the MEM at LO and NLO is its application to 
the measurement of the mass of the Z boson at the 7 TeV LHC. To this end we generate 
O(5000) events using Pythia that satisfy the following lab frame requirements, 

Pt > 15 GeV , \ru\ < 2.5 , 80 GeV < m t+e - < 100 GeV. (4.5) 

We use Pythia since it is a completely independent code to MCFM and as such is also 
independent of our new method for generating the NLO weights. In addition, Pythia output 
is at the particle level, including shower, hadronisation and underlying event models. We 
note that in Pythia we have turned off both the mass of the leptons and QED radiation, 



- 15 - 



truth 




91.05 91.10 91.15 91.20 91.25 91.30 

m z [GeV] 

Figure 4: Log-likelihoods obtained by a MEM analysis at LO (black) and NLO (red) for the 
measurement of mz at the LHC using Pythia data. Errors represent MC integration uncertainty. 

both of which ensure our transfer function assumptions remain valid. In Fig. |] we present 
the likelihoods as a function of mz for the completely inclusive case (i.e. the full data set). 
As expected we observe a parabolic function around the best fit mass. Error bars represent 
the Monte Carlo integration uncertainty and statistical uncertainties can be inferred by 
using Eq. ( |2.4|) . We observe that the truth value (mz = 91.1876 GeV) easily lies within 
the l-o" band of our best fit values, 

LO: mz = 91.170 ± 0.025 GeV NLO: mz = 91.174 ± 0.025 GeV . (4.6) 

The power of the MEM is also illustrated here, since with a data set of 0(0.1) fb _1 we are 
able to perform a measurement of the Z mass to within 25 MeV (modulo transfer function 
uncertainties). It is not surprising that the NLO and LO results are very close to one 
another since we have already observed that, for this process, the NLO and LO kinematics 
are very similar in the MEM frame. 

The results presented in Fig. || are for the full sample that includes events in which there 
is a significant amount of showered radiation. Since there is no model of this additional 
radiation in the LO MEM, one may worry that the measured value of mz depends on the 
amount of this additional radiation. We therefore present the results of a study of this 
effect in Fig. [5L where we have performed the mass measurement for a variety of cuts on 
the transverse momentum of the dilepton (Z) system, p¥f. By varying the maximum value 
of this quantity for events in our sample, we are limiting the amount of additional radiation 
(i.e. showering) present in the event. Since this veto represents an additional cut on the 
data, the size of the data sample shrinks as the maximum p^f is reduced. For this reason the 
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Figure 5: Reconstructed Z mass as a function of the upper bound on the transverse momentum 
of the dilepton system, pif. Errors represent the la deviation from the central value. Note that 
both LO and NLO calculations are performed at the same values of the cut, pip. In the plot the 
NLO points have been moved slightly to the right for clarity. 



statistical uncertainty increases at low pip, as is apparent from the uncertainties shown in 
the figure. For this observable it is clear that both the dependence on the boost and on the 
higher order corrections is small. The relative independence of the results from the amount 
of shower radiation allowed in the events illustrates that the boost method has worked well 
for this observable. This is encouraging but should not be taken as a general rule for all 
observables. The boost changes the parton fractions x a and and thus observables that 
are sensitive to such changes will become dependent on the amount of additional radiation 
in the event. In cases where imposing a jet veto is desirable, the boost (in)dependence 
should be checked by performing the measurement with a desired veto, and recalculating 
the observable with a tighter veto upon the same data set. If the two results agree within 
statistical errors then one is reassured that the shower is playing a minimal role. One may 
expect that, given its improved modeling of additional radiation, that the NLO results will 
be less sensitive to the additional radiation. 

5. The Higgs Boson search in the channel H —s- ZZ — > 4£ 

A convenient example in which to test our MEM implementation at LO and NLO is 
the Higgs search at the LHC. One of the cleanest search channels is the process H —> 
ZZ — > 4£ [39, 40] since the final state can be fully reconstructed in the detector and the 
SM backgrounds are small. With full control of the final state, with no sizeable missing 
transverse momentum or jet activity expected, this channel is a natural candidate for a 
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MEM approach. The use of the MEM in this channel has been studied in some detail in 
ref. [41], with the usual caveat of the leading order limitation. Since the NLO corrections to 
this process are large it is interesting to determine whether the MEM at NLO can improve 
upon the LO analysis. 

5.1 Event selection 

In the following examples we will select events that contain four leptons satisfying the 
following requirements, 

p^ /2 > 20 GeV , p^ M > 5 GeV , |%| < 2.0 , 

15 GeV < m a < 115 GeV , 75 GeV < m vt < 115 GeV , (5.1) 

where leptons are labelled in order of decreasing transverse momentum from i\ to £4. That 
is, we require one pair of oppositely-charged leptons to have an invariant mass within 
approximately 15 GeV of the Z mass while the invariant mass of the other pair is less 
constrained. In experimental searches the analysis cuts are typically tailored to the putative 
Higgs mass in order to better discriminate against the relevant backgrounds. However, for 
simplicity, in our studies we do not optimise the cuts in this way. Therefore the limits and 
uncertainty ranges quoted here should be taken only as a rough estimate of what can be 
achieved in a true experimental analysis. Instead, we are more interested in assessing the 
performance of the MEM at LO and NLO for a given set of cuts. 

We perform our calculation for the LHC operating at y/s = 7 TeV, with pn = fip = tuh 
in the calculation of the Higgs signal and [Ir = \ip = 2mz for the ZZ background. We 
have used the MSTW2008 PDF set [42] matched to the appropriate order in perturbation 
theory. Our NLO calculation includes the contributions from gg — > ZZ for rif = 5 massless 
flavours, using results taken from MCFM [32]. Although the interference between SM 
production of WW pairs and the Higgs signal may be phenomenologically relevant [43], in 
the ZZ — > Ai channel the corresponding interference is not expected to be important for 
a light Higgs boson since the final state is fully reconstructed. Although the interference 
effects may become non-negligible for Higgs masses above a few hundred GeV, we do not 
include such effects here. 

5.2 Higgs mass limits in the absence of the signal process 

We begin by studying the scenario in which there is no Higgs boson and the only source 
of four lepton events is pp — > ZZ — > Al production, i.e. neglecting any other source of 
backgrounds. We therefore expect to set limits on the SM Higgs, which we can do by 
defining one- and two-sigma confidence limits on its mass. As an illustrative example we 
choose to keep our definition of these contours from the Gaussian definition given in the 
previous section. Since the main thrust of this work is the study of the differences between 
LO and NLO MEM analyses, the exact definition of these limits is not crucial to our 
conclusions. 

We study the MEM using event samples that are generated at NLO parton level. 
Hadronic production of a Higgs boson through the gluon fusion mechanism is known to 



- 18 - 



have large K-factors associated with higher order corrections. In this instance this is also 
true of the continuum background. The size of the gluon flux at the LHC also results in the 
gg-initiated contributions contributing at the order of 10% of the total. Since MCFM [32] 
includes all these effects and is known to model the shape of the four lepton invariant 
mass distribution reasonably well [39], we expect that a NLO parton level simulation 
should give a reasonably realistic description of the most important features of actual four 
lepton events. We use the procedure outlined at the end of section || to directly generate 
NLO unweighted events in the MEM frame, where each of these events has the kinematic 
structure of a Born phase space point. We note that some of these events possess leptons 
with, for instance, p^ EM ' i4 < 5 GeV. Since, at LO, p^ EM = p^ b these events cannot 
pass the fiducial cuts in the LO analysis and as such are not included in the calculation 
of the likelihood. However, at NLO, the transverse momentum is not identical in the two 
frames, p^ ^ p l ^ h '. Therefore a value of p^ EM < 5 GeV can correspond to a real 
radiation contribution with p^ b > 5 GeV. As a result such events are included in the NLO 
likelihood calculation. Therefore there can be a different number of events in the LO and 
NLO data samples. This is a reflection of the fact that the NLO calculation exhibits a 
richer kinematical structure than the LO one. 

We generate 253 NLO background events, pp — > ZZ — > 4£, that pass the cuts given in 



Eq. (5JJ and then perform a MEM analysis at both LO and NLO. Due to the issue discussed 
above, only 250 of these enter the LO calculation. Our results are shown in Figure |(| where 
we present results as log-likelihood differences, log(Le/ Ls+b)- The likelihoods are for a 
signal plus background hypothesis for a given Higgs mass (Ls+b), an d for a hypothesis in 
which there is background only (Lb)- If a Higgs boson at a given ran is more probable 
than the background only hypothesis then the resulting difference is negative. However if 
the Higgs signal is less favourable than the background only hypothesis then this quantity 
is positive. Therefore a larger difference with respect to zero results in a stronger signal 
or more stringent limits. Since our sample was generated without a Higgs boson signal 
(and no chance fluctuations mimicking a signal) we observe that the difference is always 
positive. We can then proceed to set limits at the 1- and 2-a levels by requiring that 
\og(Ls / Ls+b) > 1/2 and log(Ls/ Ls+b) > 2 respectively. With these definitions we find 
the following 95% confidence exclusion range at LO, 

LO MEM exclusion: 120 GeV < m H < 380 GeV , (5.2) 

while the NLO MEM provides the more stringent exclusion limit, 

NLO MEM exclusion: 100 GeV < m H < 430 GeV . (5.3) 

Although these limits are an interesting example of the MEM at LO and NLO, we remind 
the reader that they should not be taken too literally. In particular we have not modified 
the selection cuts as a function of m#, for instance to use more efficient cuts in the heavy 
Higgs region that better maximise the signal/background ratio. In that scenario one would 
be better able to discriminate between signal and background and the presence (or lack) of 
Higgs signal events would result in stronger limits (or greater significance of a discovery). 
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Figure 6: The log-likelihood difference for background only and signal plus background, for a 
Higgs boson search in the channel, H — » ZZ* — > 4 lcptons. Positive values of the difference indicate 
that the background-only hypothesis is more likely than the signal plus background one. The blue 
and magenta lines represent the 1- and 2-tr limits respectively. 



In section 5J3 we will present a more detailed study, involving the generation of many 
pseudo-experiments, in order to investigate the true potential for setting limits on the 
Higgs boson mass by using the MEM at LO and NLO. We also note that in this section 
we are using NLO data and therefore the fact that NLO gives stronger limits is hardly 
surprising. However under the assumption that the data is better modelled by a NLO 
prediction than a LO one, we can infer that the NLO MEM might give better results than 
its LO counterpart. 



5.3 Results in the presence of a Higgs boson with ran = 125 GeV 

We now consider the case in which the SM Higgs exists and has a mass ran = 125 GeV. 
We use the same background events and cuts as in the previous section and include in 
addition three Higgs signal events that are generated from an unweighted NLO sample. 

For this mass the SM Higgs is a very narrow resonance, Tn = 0.417 x 10 -2 GeV, 
and as a result the MEM might be expected to measure the mass very precisely even 
with only a few events. This is indeed what we observe in our results, shown in Fig. [7|. 
The curves, both at LO and NLO, exhibit a very sharp valley at ran = 125 because each 
unweighted signal event has almost exactly the same mass due to the tiny width. A handful 
of events clustered at a single mass, in this case ran = 125 GeV, is clearly a very significant 
feature. In this example one signal event has p^ EM,ii = 3.8 GeV and, as already discussed, 
therefore fails the cuts in the LO MEM analysis. However, the remaining 3 leptons in this 
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Figure 7: Log likelihoods for a light Higgs boson decaying into four leptons, where we have injected 
a signal at run = 125 GeV. The dashed lines indicate the results obtained when including all 3 
signal events at NLO. The dotted line represents the NLO likelihood when the event that fails the 
cuts in the LO analysis is omitted. 



event have larger transverse momenta, thus allowing real radiation contributions to provide 
a contribution to the NLO likelihood. This enhances the differences between LO and NLO 
since the overall sample size is so small. In order to quantify the impact of the event that 
fails the cuts at LO, we present the NLO results obtained when omitting this event as 
a dotted curve in Fig. |7[ We observe that the peak value is reduced by about a third, 
indicating that although this event fails the LO cuts, it still possesses a non- negligible 
weight that the NLO calculation is able to utilise. The remaining large difference between 
LO and NLO can be understood by considering the physics in the region around my ~ 125 
GeV. As a result of our cuts there are very few background events, which means that the 
log-likelihood difference is more sensitive to the large signal l\-factor than in regions in 
which the background is sizeable. In other words, since the i'f-factor for the background 
is also large, regions in which background events occur frequently tend to spoil the NLO 
enhancement of the signal in the likelihood difference. 

In both LO and NLO cases, the MEM can measure the mass much better than the 
expected experimental resolution from the transfer functions. That is, we have stretched 
our assumption that W(x, y) = <5(x — y) well beyond its expected validity. Clearly, a 
trustworthy assessment of the accuracy of a MEM element approach in this scenario must 
be performed using a more sophisticated experimental analysis, which is beyond the scope 
of this work. We note that if the transfer functions can be determined sufficiently accurately 
(say from Drell-Yan events) then the MEM will most likely provide the best measurement 
of the Higgs mass in this channel. In addition the MEM should be very useful in the future 
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study of the light Higgs, i.e. by probing its couplings and spin. We leave such a study to 
future work. 

5.4 Results in the presence of a Higgs boson with = 425 GeV 

Finally we consider the case of a heavier Higgs boson (m# = 425 GeV), that is consequently 
much broader, Th = 37.6 GeV. Although not currently favoured experimentally this is still 
an interesting example for two reasons. Firstly the large width, compared to the case of a 
light Higgs just discussed, illustrates the behaviour of the MEM in the presence of a broad 
resonance. Secondly, although the current limits rule out a SM Higgs at approximately 
the 2(j level, they do not rule out the possibility that some other (weaker) broad resonance 
could be found in the four lepton channel. Since in this example the Higgs boson has a 
large width we expect that the likelihood results will be different to the light Higgs example 
discussed previously. 

We use the same background sample as in the previous examples and introduce four 
Higgs signal events that have been generated from an unweighted NLO sample. The 
results of a MEM likelihood analysis are presented in Fig. ^. We fit a parabola in the 
region 410 GeV < tuh < 450 GeV in order to extract a best fit mass. We observe that in 
a similar fashion to the mass measurement examples that we have already presented, LO 
and NLO analyses return similar best fit masses and uncertainties: 

mff(LO) = 428 ± 14 GeV , mf(NLO) = 427 ± 14 GeV . (5.4) 

However, since the NLO model fits the data better than the LO model the NLO result is 
more significant compared to the background-only hypothesis. 

A full understanding of the potential of the MEM for measuring heavy resonances that 
decay into leptons requires a detailed study involving many pseudo-experiments. Indeed 
since the width is so large one may often encounter individual experiments that include 
a small number of events with very large invariant mass. We expect this will result in a 
large spread of possible errors in the mass measurements. One can only perform a true 
comparison of the discriminating power of LO and NLO MEM analyses by investigating 
results obtained using an ensemble of pseudo-experiments. It is precisely this issue that 
we will investigate in more detail in the next section. 

5.5 A study with multiple pseudo experiments 

In order to investigate the differences between the MEM at LO and NLO in a more sys- 
tematic manner we repeat our analysis with multiple pseudo-experiments. We generate 
background-only data samples and study Higgs exclusion for a given Higgs mass hypoth- 
esis. We choose to investigate the case run = 200 GeV since at this value the Higgs cross 
section is near its maximum in the four lepton channel. This means that, even with a 
smaller data set than in the previous examples, we can extract useful information from an 
individual pseudo-experiment. 

Our setup is as follows. We generate 960 pseudo-experiments corresponding to un- 
weighted NLO events using the tools described in the previous sections. We attempt to 



- 22 - 



350 375 400 425 450 475 

m H [GeV] 

Figure 8: Log likelihoods for a heavy Higgs boson decaying into four leptons, where we have 
injected a signal at run = 425 GeV. 

generate the same number of background-only four lepton events as observed in the 5 fb _1 
data set of each LHC experiment. This corresponds to approximately 70 events. Any 
given pseudo-experiment will contain a random fluctuation around this central value as 
a result of the unweighting procedure. When we generate our background-only samples 
each pseudo-experiment has equal weight, regardless of the number of events entering the 
sample. Clearly a pseudo-experiment with more events will be able to set better limits and 
we have not specifically accounted for this effect. However, since we are mostly interested 
in comparing the limits obtained at LO and NLO, this caveat should not greatly affect our 
conclusions. 

In Fig. ^ we show the distribution of A = log (Lb/Ls+b) for the ensemble of pseudo- 
experiments. Since the signal at 200 GeV is relatively strong, a typical pseudo-experiment 
- that contains only background events - is able to exclude this hypothesis effectively, i.e. 
A > 0. We note that, as expected, the NLO MEM typically sets a much stronger exclusion 
than at LO (the peak in the LO distribution is in the region A ~ 5, whilst the NLO peak 
is at A ~ 9). This is also clear from the summary of exclusion limits that could be set, 
presented in Table. [|. 

Although it is interesting to examine the differences between the two orders in pertur- 
bation theory, we stress that it is not surprising that NLO is able to set stronger limits, 
since the event samples in each pseudo-experiment are generated by a NLO calculation. 
One might expect that if we had used LO unweighted samples instead, that LO would have 
set stronger limits for this Higgs mass. However, since NLO predictions typically provide 
a better description of experimental data than LO ones, there is some justification to the 
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Table 1: Percentage of pseudo experiments which set limits given by an upper bound on the 
confidence level. 



claim that the NLO MEM could set stronger limits in such a Higgs search. In addition it is 
interesting to note that the likelihoods obtained at the two orders in perturbation theory 
are significantly different, regardless of which analysis sets the better limit. Given the large 
if-factors associated with these processes, this is not unexpected. 

6. Conclusions 

The matrix element method is an analysis technique that can be used to determine param- 
eters of an underlying physics model by using a set of events that are measured experimen- 
tally. The probability that a single event in the set is described by a given model hypothesis 
can be computed from a calculation of the scattering probability within that model. Up 
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until now, the use of this technique had been limited to scattering probabilities computed 
at the leading order in perturbation theory, corresponding to Born matrix elements. 

Even at leading order, a key issue that must be addressed is the means by which a 
generic experimental event is mapped to a scattering probability. In particular, such events 
typically contain additional hadronic activity that cannot be modelled by the simplest Born 
matrix elements. In this paper we have introduced a procedure for handling this mapping in 
a consistent manner. One can combine all of the event that is not part of the desired Born 
final state into one four vector, X and then boost into a frame in which X is at rest in the 
transverse plane. This feature does not uniquely define the boost and, although the matrix 
element is a Lorentz scalar, the convolution with the parton distribution functions depends 
on the specific nature of the boost. Therefore a theoretically well-defined procedure is only 
obtained by integrating over all allowed boosts. Once this has been done we can produce 
a well defined LO weight that can be associated with each experimental event. 

We have subsequently illustrated how one can extend this method to next-to-leading 
order. The incorporation of some elements of the calculation, such as the virtual diagrams, 
is relatively straightforward since the diagrams share the same phase space as the Born 
calculation. The inclusion of the real radiation contribution is more complex and is per- 
formed by using a forward branching phase space generator. This allows one to maintain 
the exact kinematics of the Born event whilst integrating out the real radiation. We have 
used a slightly modified version of the usual Catani-Seymour dipole subtraction procedure 
in order to ensure event-by-event subtraction scheme independence. Using this generator 
we are able to define a map between all NLO events and Born phase space points. The 
final result is a method for generating a full NLO weight from a given Born phase space 
point. We note that there are some subtleties in this method that require particular care. 
For example, the difference between the lab- and MEM-frame transverse momentum can 
mean that events that are within the fiducial region in the lab frame cannot ultimately be 
included in a LO MEM analysis. At NLO, such events can be accommodated in a NLO 
MEM approach since they are accounted for by the presence of real radiation. 

We have tested the method by producing NLO likelihoods for events that contain 
electroweak final states. As a first example we considered production of lepton pairs and 
showed that one could correctly measure the mass of the Z boson using events generated 
with Pythia. In this instance we observed that the MEM frame kinematics are very similar 
at LO and NLO. For this reason we only observed small differences between the MEM at 
LO and NLO for this process. We then considered the search for the Higgs boson in the 
channel gg — > H — > ZZ* — > 4 leptons. We showed that in this case the differences between 
the LO and NLO MEM analyses could be large. By using the NLO procedure one could 
obtain better limits on Higgs exclusion and observe stronger signals should a Higgs boson 
exist. In this case our input data was generated by a NLO parton level calculation, so it 
is perhaps not surprising that the NLO MEM analysis performed better than the LO one. 
Statements regarding possible improvements in results from the MEM at NLO are therefore 
predicated on the assumption that, in general, NLO predictions describe experimental data 
better than LO ones. 

Future applications of the NLO MEM are widespread. One obvious example is the 
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measurement of the top quark mass. In addition the MEM is very useful when data samples 
are limited by statistics. Examples of measurements that would fall into this category 
include the measurement of the properties of the Higgs boson and limits on anomalous 
gauge boson couplings. We hope to extend our method to include more complicated final 
states, such as ones containing neutrinos and jets, shortly. The examples presented here 
have been implemented in a Fortran code that may be obtained from the authors on 
request. 2 
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A. An Initial State Forward Branching Phase Space Generator 

In this appendix we will discuss the generation of the forward branching phase space 
(FBPS) used in our method, additional details for which can be found in refs. [38, 44]. To 
start the derivation of the initial state FBPS generator we recall the phase space for the 
production of a heavy state Q from two initial partons, p a and p b , as, 



1 -d<S>[ D] (p a +p b 



Q) 



2vr 



(A.l) 



2s ab 2s a t 

Here we have maintained the notation used in sec. ^, where hatted momenta indicate 
the underlying Born topology and unhatted momenta represent the phase space with one 
additional parton. The D-dimensional phase space for the emission of one extra initial 
state parton with momentum p r is given by [44] , 
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We need the FBPS generator in four-dimensions 
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(A.3) 



where <f> is the azimuthal angle around the £-clxis, S X y cUld t X y are invariants defined through, 
s X y = {Px + Py) 2 and t xy = (p x —Py) 2 - While the above formula gives the phase space 
integrator, we need to derive both the integration boundaries and the explicit construction 



2 Please contact ciaran@fnal.gov for a copy of the code. 
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^ F bps = (^) i dt rb r i d<p . (a.5) 

J Sab V 



of the generated four- vectors (p a , p h and p r ) that are used in a numerical Monte Carlo 
integrator. The phase space generator starts using the input momenta p a and pb, 

p a = E a (1,0, 0,-1) , p 6 = £ 6 (l,0,0,l) , (A.4) 

with s afc = 2p a -p 6 = AE a E b . We can eliminate i ar in favour of s ab by inserting the operator 
J d s ab S(s ab + t ar + t rb — Sab) = 1 to perform the i ar integration, 

V(2tt)V y_ tm . n ^ Js ab \ S ab / JO 

The integration limits on s a b can be understood from the momentum conserving delta 
function and the requirement that t ar , t rb < 0. We will define i min shortly. Our task is 
then to construct the new momenta p a , p b and p r from the MC integration variables and 
determine the integration boundary t min . We relate s ab and t r b to our MC integration 
variable using logarithmic sampling, 

f ^ = log (±) f dr- s ab {r) = s ab s^ (A.6) 

Js ab s ab \S a b/ Jo 

/0 I — t BO ft f0 
dtbr = dti„, + / dtbr 

J — *soft 

* SOft ft ■ \ f 1 1- 

dt 6r = log I / dr t r6 (j"); i r6 (r) = - (t min ) r (t soft ) r . (A.7) 

-tmin \C B oft/ JO 

Our phase space measure is now written in terms of MC integration variables and our 
final task is to determine p a , pb and p r for use in the matrix element in terms of our new 
variables. We wish to branch one of our initial state momenta and in this example we 
choose to branch p b . In order to do so we have to give it a virtuality t rb , which we can do 
by boosting p a , 

Pa = (1 + P) Pa , Pb=Pb-f3Pa, (A. 8) 

with f3 = —trbfsab- Note that p a +Pb = Pa+Pb = Pa — Pr + Pb- This means we have added 
to the phase space generator a factor f d(3 5{(3 + t r bfs a b) that does not change the phase 
space weight. We define p a = p a = (1 + f3)p a and parametrize pb as follows, 

Pb = zE b (l,cos#, sin#cos</>, sin sin 0) , (A.9) 

where 9 is the polar angle with respect to momentum p b . Momentum conservation now 
fixes p r = Pb — Pb = Pb — Pb + PPa- To express cos# and z in terms of the integration 
variables and the input energies we calculate the invariants 

(A.10) 
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We can invert the equations to obtain. 

?2 



z 



cos 6 

By choosing 
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£ min = min(s ab - s ab , s ab (y/s - E a )/E a ) , 



we fulfill the requirement — 1 < cos 6 < 1. Remaining constraints (such as a jet-veto on p r ) 
are imposed through event vetos. 

We have thus illustrated how we have implemented the FBPS to perform both the 
integration over emitted partons and the phase space generation of p a , p b and p r for use 
in the matrix element. The input for the generator is just the Born kinematics, i.e. p a , p b 
and Q. 

B. Subtraction terms in the MEM frame 

In this appendix we discuss the modifications to the Catani-Seymour dipoles [37] needed to 
correctly ensure a one-to-one map between the integrated and unintegrated subtractions on 
an event by event basis. In this paper we consider processes with electroweak final states, 
and as such only need initial-initial dipoles. In the standard approach one would perform 
a transformation such that the emitter and spectator are kept along the beam axis, with a 
Lorentz transformation on the remaining final state particles performed in order to ensure 
momentum conservation. In our case it is essential to keep the final state particles fixed 
and instead change the momenta of the initial state partons. 

The standard Catani Seymour dipole keeps the momentum of the spectator initial 
particle b fixed, while the emitter a is rescaled by an amount x a . r , 

Par = %a,r Pa j 

Sab ~\~ S a r ~\~ S r b /t-, -, \ 

X r ,ab = • (B.l) 

Sab 

Here we have kept the same notation as the previous section, with r, a and b representing 
the emitted parton, initial state emitter and initial state spectator respectively. Hatted 
momenta still represent the underlying Born phase space - with unhatted momenta indi- 
cating the real phase space point - and in addition p now represents the dipole phase space 
point. The transformation above is given by Eqs. (5.137) and (5.138) in Ref. [37] using our 
momentum definitions. In order to ensure that p is a correct phase space point one must 
perform a Lorentz transformation (Eqs. (5.139) - (5.144) in Ref. [37]) to ensure momentum 
conservation. 

The above transformation is not ideal for our setup. This is because the Lorentz 
transformation will naturally change the underlying Born phase space point. This means 
that there will not be a one-to-one correspondence between real and virtual events and only 
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the sum over all virtual and real contributions will be well-defined. In order to maintain 
our exact map to the Born phase space p a + p& — I Q we replace Eq. ( B.l| ) by the following 
transformation, 

Par — %a,rPa ; (-^ - ^) 

Sab i S ar + S r f) . , 

X r ,ab = • (B.3) 

Sab 

Note that the transformation acts on p a , the initial state momentum of the Born phase 
space. We note that this transformation preserves momentum conservation in the trans- 
verse plane, but not in the longitudinal plane. Therefore the correct dipole phase space 
point is at a different x a and x than the original Born phase space point. Since we in- 
tegrate over these variables this is sufficient to obtain the exact mapping between virtual 
and real contributions on am event by event basis. Using our new transformation we can 
implement the usual Catani-Seymour dipole formulae (Eqs. (5.145) - (5.156) in Ref. [37]). 
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